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ABSTRACT 



We present an updated all-particle energy spectrum of primary cosmic rays in a wide range 
from 10 14 eV to 10 17 eV using 5.5 x 10 7 events collected in the period from 2000 November through 
2004 October by the Tibet-Ill air-shower array located at 4300 m above sea level (atmospheric 
depth of 606 g/cm 2 ). The size spectrum exhibits a sharp knee at a corresponding primary energy 
around 4 PeV. This work uses increased statistics and new simulation calculations for the analysis. 
We performed extensive Monte Carlo calculations and discuss the model dependences involved 
in the final result assuming interaction models of QGSJETOlc and SIBYLL2.1 and primary 
composition models of heavy dominant (HD) and proton dominant (PD) ones. Pure proton 
and pure iron primary models are also examined as extreme cases. The detector simulation was 
also made to improve the accuracy of determining the size of the air showers and the energy of 
the primary particle. We confirmed that the all-particle energy spectra obtained under various 
plausible model parameters are not significantly different from each other as expected from the 
characteristics of the experiment at the high altitude, where the air showers of the primary 
energy around the knee reaches near maximum development and their features are dominated 
by electromagnetic components leading to the weak dependence on the interaction model or the 
primary mass. This is the highest-statistical and the best systematics-controlled measurement 
covering the widest energy range around the knee energy region. 

Subject headings: cosmic rays — methods: data analysis — stars: supernovae : general 



1. Introduction 

Although almost 100 years have passed over 
since the discovery of cosmic rays, their source and 
acceleration mechanism are still not fully under- 
stood. The energy spectrum and chemical com- 
position of cosmic rays can be key information for 
probing their origin, acceleration mechanism and 
propagation mechanism. The cosmic-ray spec- 
trum has been observed by many ground-based 
experiments to resemble two power laws, having 
a form dj/dE oc E^ 1 , with 7 = 2.7 below the 
energy around 4 x 10 15 eV, and then steepen- 
ing to 7 = 3.1 above this energy (Horandel 2003). 
The change of the power index at this energy is 
called the spectral "knee" . Although the existence 
of the knee has been well established experimen- 
tally, there are still controversial arguments on 
its origin. Proposals for its origin range from as- 
trophysical scenarios like the change of accelera- 
tion mechanisms (Bcrczhko & Ksenofontov 1999; 
Stanev et al. 1993; Kobayakawa et al. 2002; Volk 2004) 
at the sources of cosmic rays (supernova rem- 
nants, pulsars, etc.), the single source assumption 
(Erlykin & Wolfendale 2005) , or effects due to the 
propagation (Ptuskin et al. 1993; Candia et al. 2002) 
inside the galaxy (diffusion, drift, escape from the 
Galaxy) to particle-physics models like the inter- 
action with relic neutrinos (Wigmans 2003) dur- 
ing transport or new processes in the atmosphere 



(Nikolsky & Romachin) during air-shower devel- 
opment. Common to all models is the prediction 
of a change of the chemical composition over the 
knee region. Direct measurements of primary cos- 
mic rays on board balloons or satellites are the 
best ways for the study of the chemical composi- 
tion, however, the energy region covered by them 
with sufficient statistics are limited to 10 14 eV. 
The energy spectrum and chemical composition of 
primary cosmic rays around the knee, therefore, 
has to be studied with ground-based air-shower 
experiments using surface array and/or detectors 
for Cerenkov light. 

Many reports have so far been made on the 
energy spectrum as well as the chemical com- 
position of primary cosmic rays. Although the 
global features of the all-particle spectrum agree 
well when we take into account the systematic er- 
rors of about 20% involved in the energy scale 
(Horandel 2003), there are still serious disagree- 
ments in the chemical composition depending on 
experimental methods, for example, the knee com- 
position obtained by the Tibet and KASCADE 
experiments can be summarized as follows. We 
have already reported the energy spectrum of 
protons and helium in the energy range from 
200 TcV to 10,000 TeV (Amenomori et al. 2000; 
Amcnomori et al. 2006a) from air-shower core ob- 
servation, suggesting a steep power index of ap- 
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proximately -3.1. This indicates that the power 
index of light component is changed from ap- 
proximately -2.7 as measured by direct observa- 
tions to -3.1 around a few hundred TeV. Hence, 
the light component should become less abundant 
at the knee and the main component responsi- 
ble to the structure of the knee must be heavier 
than helium. Furthermore, the spectral shape of 
light component seems to keep the power law in- 
stead of the exponential cutoff. On the contrary, 
the KASCADE using electron-muon size analysis 
(Antoni et al. 2005) claims that the knee in the 
all-particle spectrum is due to the steepening of 
the spectra of light elements with exponential type 
cutoff. 

The accurate measurement of the all-particle 
energy spectrum around the knee is essential to 
establish the chemical composition at this energy 
range. There is no precise measurement of the 
chemical composition around the knee region yet 
and it is in fact impossible to discriminate indi- 
vidual elements clearly by indirect observations. 
Therefore, most of the works made so far simply 
discussed the average mass < In A >. Another ap- 
proach is the unfolding of the all-particle spectrum 
using shower characteristics like electron-muon ra- 
tio, depth of the shower maximum and so on. In 
these methods, the detailed information of the all- 
particle spectrum plays an important role in de- 
termining the chemical composition. It is also 
expected that the specific features of each com- 
ponent like cutoff energy or source characteristics 
should be reflected in the shape of the all-particle 
spectrum as discussed in the single source model 
(Erlykin & Wolfendale 2005). The important fea- 
tures of the all-particle spectrum are the absolute 
intensity, the position of the knee, the difference of 
the power index before and after the knee and the 
sharpness in the size spectrum, which are deeply 
connected with the acceleration mechanism and 
the source of cosmic rays. 

The merit of the air-shower experiment in Tibet 
is that the atmospheric depth of the experimen- 
tal site (4300 m a.s.L, 606 g/cm 2 ) is close to the 
maximum development of the air showers with en- 
ergies around the knee almost independent of the 
masses of primary cosmic rays as demonstrated 
in Fig. 1 for vertically incident cosmic rays. It 
should be also noted that the number of shower 
particles is dominated by electromagnetic compo- 



nent with minor contribution of muons, whose in- 
teraction model dependence is known to be rather 
large among current interaction models leading to 
a large systematic error in the experiments carried 
out at the sea level because of the large contribu- 
tion of muons. In other words, the air shower ob- 
servation at high altitude is sensitive to the most 
forward region of the hadronic interactions in the 
center of momentum system (CMS) where high 
energy secondaries are produced, and the electro- 
magnetic component as a decay product of neutral 
pions dominate the number of shower particles, 
while it is insensitive to the central region of the 
CMS where large number of muons as decay prod- 
uct of charged pions are produced. The differences 
among current interaction models are mainly re- 
lated to the central region as seen in the prob- 
lem of the electron-muon correlation. Hence, the 
air-shower experiment in Tibet can determine the 
primary cosmic-ray energy much less dependently 
upon the chemical composition and the interaction 
model than experiments at the sea level. 

We have already reported the first result on 
the all-particle spectrum around the knee region 
based on data from 2000 November to 2001 Oc- 
tober observed by the Tibet-Ill air-shower ar- 
ray (Amenomori et al. 2003a). In this paper, we 
present an updated all-particle energy spectrum 
using data set collected in the period from 2000 
November through 2004 October. The updates are 
due to (1) increased statistics by approximately 
2.6 times, (2) use of new simulation codes and (3) 
improvement of the lateral structure function used 
for the size estimation of air showers. The previous 




Atmospheric Depth ( g/cm"2 ) 

Fig. 1. — Average transition curves of air-shower 
size induced by protons and iron nuclei for a ver- 
tical incidence. 
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result was obtained using almost the same analysis 
as used in Tibet I (Amenomori et al. 1996) except 
for the parameters which depend on the detector 
configuration. In the present paper, the simulation 
code Cosmos is replaced by Corsika with interac- 
tion models QGSJETOlc and SIBYLL2.1, which is 
now widely used in many analyses by other works 
and enables the comparison of this work with oth- 
ers easier. The third update on the structure func- 
tion is made to cover wider energy range than be- 
fore (see section 4.1.3). 

Thus, we obtained the all-particle energy spec- 
trum of cosmic rays in a wide range over 3 decades 
between 10 14 eV and 10 17 eV and the updated re- 
sult is compared with previous ones. 



Tibet III Air Shower Array (2003) 



2. Tibet experiment 

The Tibet air-shower experiment has been op- 
erated at Yangbajing (E90°31\ N30°06'; 4300 m 
above sea level) in Tibet, China, since 1990. The 
Tibet air-shower array is designed not only for 
observation of air showers of nuclear-component 
origin but also for that of high energy celestial 
gamma rays. Because of such multiple purposes, 
the detector is constructed to cover a wide dy- 
namic range for particle density covering 0.1 to 
5000 and a good angular resolution for the arrival 
direction of air showers with energy in excess of a 
few TeV being better than 1 degree. 

The Tibet-I surface array was constructed in 
1990 (Amenomori et al. 1992) using 65 plastic 
scintillation detectors placed on a lattice with 15 
m spacing. This array was gradually expanded 
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Fig. 2. — Schematic view of the FT w/D-detector. 



i i □ FT Detector D-PMT(249) 

,5m ;:■ Density Detector ( 28} 

Fig. 3. — Schematic view of the Tibet-Ill ar- 
ray operating at Yangbajing. The Tibet-Ill array 
consists of 761 FT detectors and 28 D detectors 
around them. In the inner 36,900 m 2 , FT detec- 
tors are deployed at 7.5 m lattice intervals, among 
761 FT counters, 249 sets of detectors are also 
equipped with D-PMT in addition to FT-PMTs. 
Open-white squares: FT detectors with FT-PMT; 
Open-black squares: FT detectors with FT-PMT 
and D-PMT; Open circles: density detectors only 
with D-PMT. 
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to the Tibet-II (1994) and Tibet-Ill (1999) array. 
At present, it consists of 761 fast timing (FT) 
counters and 28 density (D) counters surrounding 
them. In the inner 36,900 m 2 , FT counters are 
deployed at 7.5 m lattice intervals. All the FT 
counters are equipped with a fast-timing photo- 
multiplier tube (FT-PMT ; Hamamatsu H1161) 
measuring up to 15 particles. Among the 761 FT 
counters, 249 sets of detectors (with interval of 
15 m) are also equipped with density photomulti- 
plicr tube (D-PMT ; Hamamatsu H3178) of wide 
dynamic range measuring up to 5000 particles in 
addition to FT-PMTs, so that UHE cosmic rays 
with energy above the knee can be observed with 
a good accuracy. 

Each counter has a plastic scintillator plate (BI- 
CRON BC-408A) of 0.5 m 2 in area and 3 cm in 
thickness. A lead plate of 0.5 cm thick is put 
on the top of each counter as shown in Fig. 2 
in order to increase the counter's sensitivity by 
converting photons in an electromagnetic shower 
into electron-positron pairs (Bloomer et al. 1988; 
Amcnomori et al. 1990). The recording of signals 
is made on time and charge information for the 
FT-PMTs, while only the charge information for 
the D-PMTs. The D counters surrounding the in- 
ner array are also equipped with both FT-PMT 
and D-PMT, where only the charge information 
of both PMTs are recorded. An event trigger sig- 
nal is issued when any fourfold coincidence occurs 
in the FT counters recording more than 0.6 parti- 
cles. Fig. 3 is the schematic view of the Tibet-Ill 
array. 

The primary energy of each event is determined 
by the shower size N e , which is calculated by 
fitting the lateral particle density distribution to 
the modified NKG structure function (see section 
4.1.3). The air-shower direction can be estimated 
with an inaccuracy smaller than 0.2° at energies 
above 10 14 eV, which is calibrated by observ- 
ing the Moon's shadow (Amenomori et al. 2003b). 
We used the data set obtained during the period 
from 2000 November through 2004 October. The 
effective live time used for the present analysis is 
805.17 days. 

3. Simulation 

Monte Carlo simulation (MC) plays an impor- 
tant role in air-shower experiments since most of 



the methods of the analysis are developed so that 
they can reproduce the inputs of simulated events 
like the primary energy, the location of the shower 
axis, the arrival direction and so on. Even the 
most basic quantity like the number of particles 
arriving to a detector should be 'defined' through 
MC because we do not measure the number of 
particles but the charge of PMT-output which is 
not simply proportional to the number of charged 
particles entering into the detector if we take into 
account the contribution of the electromagnetic 
processes by photons inside the detectors. An- 
other example of the role of MC is to define the 
effective area of the shower array which should be 
determined to avoid the erroneous counting of the 
events whose shower axises are dropping outside of 
the effective area. Therefore detailed MC calcula- 
tions are needed on air-shower generation in the 
atmosphere and on the detector response. Con- 
sequently, the final result inevitably depends on 
the interaction model and on the primary compo- 
sition model in MC. This is the main source of the 
systematic errors involved in the air-shower exper- 
iment and we try to show them explicitly in the 
present work. 

A full Monte Carlo (MC) simulation has been 
carried out on the development of air showers in 
the atmosphere and also on the detector response 
of the Tibet-Ill array. The simulation code COR- 
SIKA (version 6.204) including QGSJETOlc and 
SIBYLL2.1 interaction models (Heck et al. 1998) 
is used to generate air-shower events. All shower 
particles in the atmosphere are traced down to the 
minimum energy of 1 MeV without using thinning 
method. 

Although the chemical composition of the pri- 
mary particles around the knee region is not well 
established yet, we have to assume it in the sim- 
ulation. The simplest way to bracket all possi- 
bilities is to assume pure proton and pure iron 
primaries. Since it is almost evident that such 
assumptions are not realistic and lead to unac- 
ceptable results showing disagreement with the 
direct observations, these results will be men- 
tioned as extreme cases. More realistic treatment 
of the chemical composition is to extrapolate the 
known composition at low energies measured by 
direct observations. The uncertainty in extrap- 
olating to the high energy range can be treated 
by bracketing the reported results on the compo- 
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Table 1: Fractions of the proton(P), helium(He), 
CNO(M), NaMgSi(H),SClAr(VH) and iron(Fc) 
components in the assumed primary cosmic- 
ray spectrum of the HD and PD models 
(Amenomori et al. 2000). 



HD model 


10 i4 -10 if> eV 


10 lb -10 ib eV 


10 ie -10 iY eV 


P 


22.6% 


11.0% 


8.1% 


He 


19.2% 


11.4% 


8.4% 


M(CNO) 


21.0% 


22.6% 


17.8% 


H(NaMgSi) 


9.0% 


9.4% 


8.1% 


VH(SClAr) 


5.6% 


6.2% 


5.8% 


Fe 


22.2% 


39.1% 


51.7% 


PD model 


10 i4 -10 lf> eV 


10 lf> -10 ib eV 


10 ib -10 iY eV 


P 


39.0% 


38.1% 


37.5% 


He 


20.4% 


19.4% 


19.1% 


M(CNO) 


15.2% 


16.1% 


16.5% 


H(NaMgSi) 


9.4% 


9.9% 


10.2% 


VH(SClAr) 


5.8% 


6.2% 


6.3% 


Fe 


9.4% 


9.9% 


10.2% 



sition study around the knee. In order to exam- 
ine the composition dependence involved in the 
all-particle spectrum, we used two kinds of the 
mixed composition models. One is based on the 
dominance of heavy component around the knee 
as reported by works (Amenomori et al. 2006a; 
Amenomori et al. 2000; Ogio 2004) and this com- 
position is called HD model. Another is based 
on the dominance of light component (p and 
He) as reported by works (Antoni et al. 2005; 
Raino 2004; Fowler 2001) and called PD model. 
The energy spectra of individual mass groups in 
HD model and PD model are shown in Fig. 4(a) 
and Fig. 4(b), respectively. Table 1 shows their 
fractional contents at given energies. In total, four 
kind of the primary composition, namely, pure 
proton, pure iron, HD model and PD model are 
used in the simulation with the minimum primary 
energy of 50 TeV. 

All secondary particles are traced until their 
energies become 1 MeV in the atmosphere. The 
shower axis was placed on Tibet array at random 
within a radius of 100 m from the center of the ar- 
ray. In order to treat the MC events in the same 
way as experimental data analysis, simulated air- 
shower events were input to the detector with the 
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Fig. 4. — Primary cosmic ray composition for 
(a) the HD model and (b) the PD model. The 
all-particle spectrum, which is a sum of each 
component, is normalized to the Tibet data, 
and they are compared with other experiments. 
PROTON satellite (Grigorov et al. 1971), 
AKENO (Nagano et al. 1984),JACEE 

(Asakimori et al. 1998), RUN JOB 
(Apanasenko et al. 2001). 
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same detector configuration as the Tibet-Ill array 
with use of Epics code (ver. 8.64) (Kasahara) to 
calculate the energy deposit of these shower parti- 
cles. Experimentally, the number of charged parti- 
cles is defined as the PMT output (charge) divided 
by that of the single particle peak, which is deter- 
mined by a probe calibration using cosmic rays, 
typically single muons. For this purpose, a small 
scintillator of 25 cm x 25 cm x 3.5 cm thick with a 
PMT (H1949) is put on the top of the each detec- 
tor during the maintenance period. This is called 
a probe detector and is used for making the trig- 
ger of the each Tibet-Ill detector. The response 
of each detector is calibrated every year through 
probe calibration. In the simulation, these events 
triggered by the probe detector was also exam- 
ined by a MC calculation, in which the primary 
particles were sampled in the energy range above 
the geomagnetic cutoff energy at Yangbajing (>10 
GeV), and all secondary particles which pass the 
probe detector and the Tibet-Ill detector simulta- 
neously were selected for the analysis. Since the 
value of PMT output is proportional to the energy 
loss of the particles passing through the scintilla- 
tor, the peak position of the energy loss distribu- 
tion corresponds to the experimental single peak 
of the probe calibration. 

According to the MC, the peak position of 
the energy loss in the scintillator is calculated 
as 6.11 MeV (the details are written in the pa- 
per (Amenomori et al. 2007)). We can then cal- 
culate the number of charged particles for each 
detector hit as the total energy loss in each scin- 
tillator divided by 6.11 MeV instead of counting 
the number of charged particles arriving to the 
detector in MC events. We confirmed that the 
shape of the energy loss distribution, which is de- 
termined by probe-calibration simulation, shows 
a reasonable agreement with the charge distribu- 
tion of the experimental data as shown in Fig. 5, 
where the proportionality between the energy loss 
AE and the PMT output charge Qi is assumed 
AE, where ki is a proportional con- 
stant depending on given detector typically being 
around 4 pC /MeV. Thus, all detector responses in- 
cluding muons and the materialization of photons 
inside the detector are taken into account. The 
total number of charged particles of each event 
(hereafter, we call this an " air shower size (Ne)") 
was estimated using the modified NKG lateral dis- 



tribution function which is tuned to reproduce 
the above defined number of particles using the 
Monte Carlo simulation under our detector con- 
figurations. The number of MC events, typically 
for QGSJET+HD, is as follows. About 10 million 
of air showers are generated with primary energy 
above 50 TeV. After imposing the selection crite- 
ria described in sec. 4.2, the surviving number of 
events is about 5 million, among which 1.9 mil- 
lion events belong to the unbiased energy region 
corresponding to E0>100 TeV. Almost the same 
number of MC events are obtained for other mod- 
els to compare with each other. 

4. Analysis 

4.1. Reconstruction of air-showers 

An example of the shower profile obtained by 
Tibet-Ill array is shown in Fig. 6(a) and Fig. 6(b) 
which represent the map of arrival time and par- 
ticle density of shower particles, respectively. Al- 
though the Tibet array has quite low energy 
threshold of a few TeV for the purpose of the 
celestial gamma ray observation, the detection ef- 
ficiency for the nuclear component including iron 
nuclei is not sufficient at low energy region. Ad- 
ditional event selection condition is required for 
the unbiased detection of all particles and for the 
capability of the lateral density fitting. The fol- 
lowing condition was applied on the selection of 
the events for the all-particle spectrum analysis. 



N D > 10 with n v > 5 



(1) 



where Njj expresses the number of detectors hit 
and n p the number of particles per a detector. 
This condition satisfies the requirements for unbi- 
ased analysis in the energy range above 100 TeV 
as described below. 

4-1.1. Determination of the core position. 

The core position of each air-shower (X core , 
Y core ) is estimated by using the following equation 



(-^ core j ^cor- 



Zpi w Zpi 

where pi is the particle density at the i-th detector 
and the weight w is an energy dependent param- 
eter varying between 0.8 and 2.0. It is confirmed 
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Fig. 5. — Charge distribution in a detector mea- 
sured by probe calibration (see the text). In or- 
der to compare the charge distribution with the 
simulation on the energy loss in a scintillator, the 
MC result is adjusted by multiplying a constant 
to meet with the same peak position as the exper- 
iment. The fluctuation of the number of photons 
in scintillation light is taken into account with the 
normal distribution in MC. 



(a) 




that the mean error of the core position can be es- 
timated as 5 m by reconstructing MC events (see 
Fig. 7). Lower energy event selection condition 
than eq. (1) leads to worse core resolution which 
makes lateral density fitting difficult. 

4-1.2. Determination of the arrival direction. 

The arrival direction of the air shower is es- 
timated using the time signal measured by the 
761 FT (fast-timing) counters. The shape of the 
shower front is assumed to be a reverse-conic type 
as shown in Fig. 8. Direction cosine of the shower 
axis is determined by using a method of least 
squares in which the difference is minimized be- 
tween the arrival time signals of each detector 
and the expected values on the assumed cone with 
given direction cosine. 

An experimental check of the angular resolu- 
tion by this method is made by the observation of 
the Moon's shadow (Amenomori et al. 2003b) us- 
ing large statistics of low energy events (>3 TeV) 
and also a so-called even-odd method, which many 
authors have been using for estimating the angular 
resolution (Amenomori et al. 1990). The recon- 
struction of the high energy MC events assures us 
that the mean error of air-shower direction can be 
estimated as 0.2° at energies above 10 14 eV (see 
Fig. 9). 

4-1.3. Estimation of the shower size. 

The lateral density distribution is corrected to 
the inclined plane perpendicular to the shower axis 
and used for the shower size estimation. In this 




Fig. 6. — An example of the map of (a) arrival 
time of shower particles, (b) particle density, ob- 
tained by Tibet-Ill array. 
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Fig. 7. — The distribution of core-position error. 
The mean error of core position can be estimated 
as 5 m. Shower selection criteria : E > 100 TeV, 
sec(i9) < 1.1 and the core position located at the 
inner 135 m x 135 m of the array. 
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Fig. 8. — Determination of air-shower direction. 
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Fig. 10. — The numerical values of a and b are 
plotted as a function of s, where original defini- 
tions of a(s) and b(s) in NKG function are shown 
by the dotted lines, and the open circles denote the 
averaged MC data using QGSJETOlc+HD model 
which are fitted by empirical formulae shown by 
red lines. See the text. 
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Fig. 9. — The distribution of opening angle be- 
tween true and estimated arrival directions. The 
mean error of air-shower direction can be esti- 
mated as 0.2°. Shower selection criteria : Eq > 
100 TeV, sec(6>) < 1.1 and the core position lo- 
cated at the inner 135 m x 135 m of the array. 



work, the determination of the lateral distribu- 
tion function of shower particles is very important, 
since the total number of charged particles in each 
event is estimated by fitting this function to the 
experimental data. Using the Monte Carlo data 
obtained under the same conditions as the exper- 
iment, we found that the following modified NKG 
function can be fitted well to the lateral distribu- 
tion of shower particles under a lead plate of 5 mm 
thickness: 

C(s) = 2TrB(a(s) + 2, -b(s) - a{s) - 2) (4) 

where r m = 30 m, and the variable s corresponds 
to the age parameter, N e is the total number of 
shower particles and B denotes the beta function. 
The original meaning of r rn in NKG formula is a 
Moliere unit, being 130 m at Tibet altitude, how- 
ever, we treat r m as a unit scale of the lateral dis- 
tribution suitable to describe the structure of the 
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air showers observed by Tibet-Ill array, whose ef- 
fective area is 135 m x 135 m. The functions a(s) 
and b(s) are determined as follows. In CORSIKA 
simulation, the shower age parameter s is calcu- 
lated at observation level by fitting the number 
of particles to a function for the one dimensional 
shower development. It may be possible to as- 
sume that air showers with the same shower age s, 
are in almost the same stage of air-shower devel- 
opment in the atmosphere, i.e. they show almost 
the same lateral distribution for shower particles 
irrespective of their primary energies. The lateral 
distribution of the particle density obtained by the 
simulation with carpet array configuration is nor- 
malized by the total number of particles which is 
derived from the total energy deposit in an in- 
finitely wide scintillator. These events are then 
classified according to the stage of air-shower de- 
velopment using the age parameter and they are 
averaged over the classified events. The fitting of 
the equation (3) to the averaged MC data is made 
to obtain the numerical values a and b. Thus, 
we can obtain the behavior of a and & as a func- 
tion of s as shown in Fig. 10(a) and Fig. 10(b), 
where original definitions of a(s) and b(s) in NKG 
function are shown by dotted lines. Although our 
result shows different dependences of a and b on 
s from the original NKG function, it is confirmed 
that the lateral distribution of the shower particles 
is better reproduced by our formula (see Fig. 11). 
This expression is valid in the range of s = 0.6 ~ 
1.6, sec9 < 1.1 and r = 5 ~ 3000 m. Two interac- 
tion models of QGSJETOlc and SIBYLL2.1 and 
four primary composition models of pure proton, 
pure iron, HD and PD are used independently to 
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Fig. 11. — Lateral density distribution of the 
charged particle obtained with use of the carpet 
simulation. The shower size Ne is better repro- 
duced by our modified NKG function. 



Table 2: The shower-size resolution are summa- 
rized for the events induced by primary parti- 
cles of E ~ 100 TeV and E Q ~ 1000 TeV with 
sec(0) < 1.1 and the core position located at in- 
ner 135 m x 135 m of the array, based on the 
QGSJET+HD, QGSJET+PD, SIBYLL+HD and 
SIBYLL+PD model. 
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determine the functions a(s) and b(s) and used in 
the analysis described below. Other details are 
described in Ref . ( Amenomori et al. 2007). 

Based on the Monte Carlo simulation, the cor- 
relation between the true shower size (true size) 
and the estimated shower size (fit size) is demon- 
strated in Fig. 12(a) and Fig. 12(b). Here, the true 
shower size means the number of particles calcu- 
lated for a carpet array while estimated shower 
size is for the real Tibet-Ill array using the mod- 
ified NKG function mentioned above. As seen in 
these two figures, a good correlation is obtained 
between the true shower size and the estimated 
shower size. The systematic deviation of less than 
1 % seen around 100 TeV of Fig. 12(a) shows 
that we need finer tuning of modified NKG func- 
tion at low energies and this error is finally cor- 
rected. The shower size is well reproduced with 
standard deviation of 5% around the primary en- 
ergy of 1000 TeV with 1.0 < sec(0) < 1.1 based 
on the QGSJET+HD model. The shower size res- 
olutions estimated are summarized for the events 
with different simulation model combinations in 
Table 2. 

4.2. Data selection 

The following event-selection criteria are adopted 
in the present analysis. 

1) More than 10 detectors should detect a signal of 
more than five particles per detector, as mentioned 
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Fig. 12. — (a) The correlation between the true 
shower size and the estimated shower size (fit size) . 
(b) The shower size resolution is estimated to be 
5% around the primary energy of 1000 TeV based 
on the QGSJET+HD model. Shower selection cri- 
teria : E > 100 TeV, sec(0) < 1.1 and the core 
position located at the inner 135 m x 135 m of 
the array. 



in eq.(l). 

2) In order to minimize the primary mass depen- 
dence on the air-shower size at Yangbajing alti- 
tude, the zenith angle (6) of the arrival direction 
of air showers should be smaller than 25°, or sec 9 
< 1.1. 

3) The rejection of the events falling outside the 
effective of the array is made by estimating 
the core position using eq. (2) with high weight 
of w = 8.0 for the particle density. Then, we im- 
posed that the core position should be inside the 
innermost 135 m x 135 m area (18225 m 2 ). This 
area is chosen with the use of MC events, so that 
the following two cases are just canceling out each 
other, namely the number of events originally in- 
side of this area but falling outside after event re- 
construction equals to the number of events in the 
opposite case. It is confirmed by simulations that 
Out-In events occupy about 10% of all selected 
events (see Fig. 13), and the energy spectrum of 
Out-In and In-Out events are almost the same, 
hence the effect of the difference between them to 
the all selected events is less than 2% at each en- 
ergy bin. 

4.3. Trigger efficiency 

Our simulations confirmed that the air show- 
ers induced by primary particles with E > 100 
TeV can be fully detected without any bias under 
the above mentioned criteria as shown in Fig. 14. 
The total effective area S x ft is then calculated 
to be 10410 m 2 -sr for all primary particles with 
Eq > 100 TeV using inner area of 135 m x 135 m 
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Fig. 13. — Out-In events occupy about 10% of 
all selected events. The energy spectrum of Out- 
In and In-Out events are almost the same, being 
the effect of the difference between them to the all 
selected events less than 2% at each energy bin. 
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and solid angle for sec(9) < 1.1. For the calcula- 
tion of the absolute intensity, the inclination effect 
due to a flat surface detector is taken into account 
by correcting the density of the observed events 
into that for a plane perpendicular to the shower 
axis. For the operation period from 2000 Novem- 
ber through 2004 October, the effective live time 
T is calculated as 805.17 days. The total number 
of air showers selected under the above conditions 
is 5.5 x 10 7 events after the inclination correction. 

5. Results and Discussions 

5.1. The model dependence of the size 
spectrum 

As a first step in checking the model depen- 
dence, the difference of the size spectra derived by 
the lateral fitting is examined using the structure 
functions based on five models (QGSJET interac- 
tion model with four primary composition models 
and SIBYLL+HD model). Shown in Fig. 15 is the 
model dependence of the air-shower size spectrum 
of nearly vertical air showers, ft is seen that the 
model dependence of the air-shower size is small 
(less than 5% in absolute intensity on the primary 
composition model, and also less than 5% on the 
hadronic interaction models). 

5.2. All-particle energy spectrum 

5.2.1. Determination of the primary energy 

In Fig. 16 we show the scatter plots of the pri- 
mary energy E , and the estimated shower size N e 
based on the QGSJET+HD model. 

Fig. 17 shows the correlations between the 
estimated shower size N e and the primary en- 
ergy E for QGSJET + 4 primary models and 
SIBYLL+HD model. In this figure, the results of 
pure composition models are also shown just to 
help the understanding of the composition depen- 
dence in the energy determination. Using these 
pure composition model will result in remarkably 
different primary energy spectrum as shown later. 
The correlation between the estimated shower size 
N e and the primary energy Eq for sec 6 < 1 . 1 can 
be well fitted using a following conversion func- 
tion for each interaction models and composition 
models. 
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Fig. 14. — The trigger efficiency of air showers. In 
the case of Njj > 10, n p > 5 and sec(0) < 1.1 with 
the core position located at the inner 135 m x 135 
m of the array, the air showers induced by protons 
or irons with E a > 100 TeV can be fully detected 
without any bias, (a) QGSJET and (b) SIBYLL. 
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Fig. 15. — The model dependence of the size spec- 
trum of nearly vertical air showers (sec(6>) < 1.1). 
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Fig. 16. — Scatter plots of the primary energy Eq 
and the estimated shower size N e based on the 
QGSJET+HD model. 



Wle + 03 



— QGSJET+lron 




- QGSJET- 


-HD 




— SIBYLL + 


HD 




- - QGSJET- 


-PD 




- - QGSJET- 


-Proton 











1e+D5 1e+06 1e+07 1e+08 

Shower Size Ne 



Fig. 17. — The correlations between the estimated 
shower size N e and the primary energy Eq for 
given models. 
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Fig. 18. — The differential energy spec- 
tra of all particles obtained by the 
present work using 5 models and they 
are compared with direct observations. 
JACEE (Asakimori et al. 1998), RUN JOB 

(Apanasenko et al. 2001), Grigorov 
(Grigorov et al. 1971) 



Eo = ax{ Tok } [TeV] ' (5) 

where the numerical values of a and b in eq.(5) are 
summarized in Table 4. These approximations are 
valid for nearly vertical air showers (sec(#) < 1.1 
) at Yangbajing altitude. 

As seen in Fig. 17, the difference of the conver- 
sion factor between two mixed composition models 
of HD and PD is not significant and almost dis- 
appears above a few times 1000 TeV in spite of 
the large difference of the fractional abundances 
of chemical components. From the comparison be- 
tween QGSJET+HD and SIBYLL+HD, one can 
see that the dependence on the interaction model 
is very small. 

It is also noted that the energy resolutions in 
different interaction models and two mixed com- 
position models are very close to each other. The 
differences are estimated as 36% and 17% at ener- 
gies around 200 TeV and 2000 TeV, respectively, in 
the case of QGSJET+HD model. Corresponding 
values for SIBYLL+HD model are 38% and 19% 
and that for QGSJET+PD model are 39% and 
19%. It is commonly understood that there is an 
energy dependent overestimation problem of the 
flux due to the steep power index of the cosmic ray 
energy spectrum when the error of the estimated 
energy increases with decreasing primary energy. 
It may be, however, worthwhile to note here that 
this effect is already included in our method of the 
energy determination using the conversion func- 
tion, because the primary energy at a given shower 
size is determined including the contribution of all 
possible primary energies which leads to a smaller 
value than the case of without fluctuation reflect- 
ing the larger population of low energy primaries 
than high energies. To avoid methodical system- 
atic error, the reproducibility of the primary flux 
was carefully examined using MC events and no 
significant deviation was found between MC input 
spectrum and the reconstructed spectrum. 

5.2.2. Energy spectrum of all-particles and the 
knee parameters 

The all-particle energy spectrum of primary 
cosmic rays in a wide range over 3 decades be- 
tween 1 x 10 14 eV and 1 x 10 17 eV is shown in 
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Fig. 19. — The differential energy spectra of all 
particles obtained by the present work using mixed 
composition models compared with other experi- 
ments. JACEE (Asakimori et al. 1998), RUN JOB 
(Apanasenko et al. 2001), Grigorov 
(Grigorov et al. 1971), KASCADE 
(Antoni et al. 2005), CASA-MIA 
(Glasmacher et al. 1999), AKENO(1984) 
(Nagano et al. 1984), Tibet-III(ICRC2003) 
(Amenomori et al. 2003a). 
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Fig. 18 for 5 models. In this figure, it is impor- 
tant to note that the position of the knee is clearly 
seen at the energy around 4 PeV irrespective of the 
model used. 

The model dependence of the conversion func- 
tion shown in Fig. 17 may lead to different shapes 
of the all-particle energy spectrum. The model 
dependence on the primary composition can be 
checked by comparing the results of QGS JET+HD 
and QGSJET+PD. Although the composition is 
fairly different between the HD and PD model at 
energies above 10 15 eV, the difference of the abso- 
lute intensity is 20% at most between the two mod- 
els and decreases with increasing primary energy. 
As mentioned above, it may be helpful to show the 
composition dependence by the pure component 
assumption. The difference between pure proton 
and pure iron primary model becomes large at the 
lower energy region of the range shown in this fig- 
ure, where the difference of the intensities between 
two models exceeds a factor of 3 at 10 14 eV, al- 
though these results do not belong to the argu- 
ment for the reality at least in the energy range 
lower than a few times 10 14 eV. It should be noted 
that the mixed composition models give more re- 
alistic extension of the direct observations and the 
composition model dependence almost disappears 
above 10 16 eV. This is the remarkable character- 
istics of the Tibet experiment. 

The interaction model dependence is seen 
by comparing the results of QGSJET+HD and 
SIBYLL+HD. This comparison shows that the 
shapes of the spectrum from both interaction 
models are almost the same and the difference in 
the absolute intensity is within 10%. In Fig. 19, 
we show the results using the mixed composition 
model in comparison with other works including 
our previous work presented in ICRC2003. Our 
results seem to lead to the data of direct obser- 
vation smoothly in the low energy side. Fig. 20 
shows the higher energy part above the knee to 
compare with the surface array experiments done 
at the highest energy region. It is very interesting 
to see how the smooth extension of our spectrum 
above 10 16 eV is connected with the data at the 
highest energy region. 

The intensities of all-particle energy spectra 
measured by the Tibet-Ill array are also posted 
in Table 5 and the summary of the measured all- 
particle energy spectrum and knee parameters are 
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Table 3: The summary of the knee parameters. 
The symbol 71 is the best fitted-index for the en- 
ergy range 100 TeV < Eo < 1 PeV, and 72 is for 
the energy above 4 PeV. 



Model 


knee position 
(PeV) 




Index of 

spectrum 


QGSJET 


(4.4 ± 0.1) 


7i 


= -2.81 ± 0.01 


+Iron 




72 


= -3.21 ± 0.01 


QGSJET 


(4.0 ± 0.1) 


7i 


= -2.67 ± 0.01 


+HD 




72 


= -3.10 ± 0.01 


QGSJET 


(3.8 ± 0.1) 


7i 


= -2.65 ± 0.01 


+PD 




72 


= -3.08 ± 0.01 


SIBYLL 


(4.0 ± 0.1) 


7i 


= -2.67 ± 0.01 


+HD 




72 


= -3.12 ± 0.01 


QGSJET 


(3.4 ± 0.1) 


7i 


= -2.60 ± 0.01 


+Proton 




72 


= -3.03 ± 0.01 



listed in Table 3, where 71 is the best fittcd-indcx 
for the energy range 100 TeV < Eq < 1 PeV, and 
72 is for the energy above 4 PeV. 

6. Summary 

We have analyzed the air shower data set col- 
lected in the period from 2000 November through 
2004 October with the Tibet-Ill air-shower array, 
using the new simulation code CORSIKA, and ob- 
tained the all-particle energy spectrum of primary 
cosmic rays in a wide energy range over 3 decades 
between 10 14 eV and 10 17 eV. The knee of the pri- 
mary spectrum is clearly observed and its position 
is located at the energy around 4 PeV. 

The advantage of the Tibet experiment at high 
altitude is first that the primary energy for the 
unbiased detection of air showers is sufficiently 
low for the purpose of the measurement around 
the knee, and second that the energy determina- 
tion is insensitive to the number of muons which 
is dependent on hadronic interaction model and 
the chemical composition. In order to quanti- 
tatively confirm these characteristics, the model 
dependence on the primary chemical composition 
was estimated in terms of the two mixed chemi- 
cal composition models of HD (heavy dominant) 
and PD (proton dominant), which are extrapo- 
lated from the direct observations with different 
scheme of the fractional contents of the individ- 
ual elements, together with the extreme cases of 
pure proton and pure iron. The interaction model 



dependence was discussed using QGSJETOlc and 
SIBYLL2.1 interaction models as they are widely 
used by other works. It was shown that the air 
showers induced by primary energies above 100 
TeV are fully detected for all kind of primary par- 
ticles. The systematic errors due to the above 
mentioned model dependences are shown to be 
within a few tens % in the energy range below the 
knee, i.e., 20% in chemical composition models be- 
tween HD and PD, and 10% in interaction models 
between QGSJETOlc and SIBYLL2.1. The model 
dependence is decreasing with increasing primary 
energy, and it almost disappears above 10 16 eV. 
Although these estimates are limited for the cho- 
sen models, other choice of more adequate models, 
if any, will not change the result of this work dras- 
tically because of the weak dependence on the 
model used. The uncertainty due to the interac- 
tion model will even decrease after the measure- 
ment of the CMS forward region by the forth- 
coming experiment LHCf (Bonechi et al. 2006; 
Sako et al. 2007). This is the highest-statistical 
and the best systematics-controlled measurement 
covering the widest energy range around the knee 
energy region. 

As discussed above, the main uncertainty left 
in the all-particle spectrum is related to the chem- 
ical composition. This work will be extended 
to the analysis of the air showers with large ar- 
rival zenith angle which can show another fea- 
tures of the air-shower development in the at- 
mosphere and provides the information about 
the chemical composition (to be published else- 
where). As mentioned before, we have already 
reported a paper on the proton and helium spec- 
tra around the knee (Amenomori et al. 2006a; 
Amenomori et al. 2000) derived from the hybrid 
experiment using the air-shower core detector, 
which is sensitive to the showers of light element 
origin like proton or helium by selecting the high 
energy core. From the observed steep power in- 
dex and the low intensities of proton and helium 
spectra, the dominance of the heavy elements 
were suggested by the hybrid experiment, how- 
ever, the statistics were limited due to the high 
threshold. In the very near future, we will start 
a new high-statistics hybrid experiment in Tibet 
(Huang et al. 2005) to clarify the main compo- 
nent of cosmic rays at the knee. The core detector 
will consist of 400 burst detectors located at the 
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center of the Tibet-Ill array in a grid with the 
detector interval of 3.75m. The burst detectors 
measure the high energy electromagnetic cascade 
of the energy above 30 GeV developed in lead plate 
of 3.5 cm thick by shower-core particles. The new 
experiment is able to observe the air-shower cores 
induced by heavy components around and beyond 
the knee, where direct measurements are inacces- 
sible because of their extremely low fluxes. The 
first observation of the iron spectrum in the knee 
region is expected in this new experiment. 
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Shower Arrays has been performed under the aus- 
pices of the Ministry of Science and Technology 
of China and the Ministry of Foreign Affairs of 
Japan. This work was supported in part by Grant- 
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from the Ministry of Education, Culture, Sports, 
Science and Technology, by Grants-in-Aid for Sci- 
ence Research from the Japan Society for the Pro- 
motion of Science in Japan, and by the Grants- 
in-Aid from the National Natural Science Foun- 
dation of China and the Chinese Academy of Sci- 
ences. The authors thank J. Kota for reading the 
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Table 4: The parameters in the conversion function 
eq. (5). 



Model 


EO < 10 ib eV 


E0 > 10 ib cV 




a 


b 


a 


b 


QGSJET+Proton 


1.195 


0.959 


1.195 


0.959 


QGSJET+HD 


1.872 


0.924 


1.348 


0.953 


QGSJET+PD 


1.583 


0.933 


1.348 


0.951 


SIBYLL +HD 


1.968 


0.913 


1.323 


0.951 


QGSJET+Iron 


3.915 


0.851 


3.915 


0.851 
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Table 5: The intensity of all-particle energy spec- 
trum measured by Tibet-Ill array based on the 
QGSJET+HD, QGSJET+PD and SIBYLL+HD 
models. 
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